#pragma once

#include <algorithm>
#include <eigen3/Eigen/Dense>
#include <functional>
#include <iostream>
namespace eg = Eigen;
const eg::Vector<double, 5> gaussp5 = {-0.90617985, -0.53846931, 0.0,
                                       0.53846931, 0.90617985};
const eg::Vector<double, 9> gaussp9 = {-0.96816024, -0.83603111, -0.61337143,
                                       -0.32425342, 0.,          0.32425342,
                                       0.61337143,  0.83603111,  0.96816024};
const eg::Vector<double, 7> gaussp7 = {-0.94910791, -0.74153119, -0.40584515,
                                       0.,          0.40584515,  0.74153119,
                                       0.94910791};
const eg::Vector<double, 11> gaussp11 = {
    -0.97822866, -0.8870626, -0.73015201, -0.51909613, -0.26954316, 0.,
    0.26954316,  0.51909613, 0.73015201,  0.8870626,   0.97822866};

const eg::Vector<double, 11> weight11 = {
    0.05566857, 0.12558037, 0.18629021, 0.23319376, 0.26280454, 0.27292509,
    0.26280454, 0.23319376, 0.18629021, 0.12558037, 0.05566857};
const eg::Vector<double, 9> weight9 = {0.08127439, 0.18064816, 0.2606107,
                                       0.31234708, 0.33023936, 0.31234708,
                                       0.2606107,  0.18064816, 0.08127439};
const eg::Vector<double, 7> weight7 = {0.12948497, 0.27970539, 0.38183005,
                                       0.41795918, 0.38183005, 0.27970539,
                                       0.12948497};
const eg::Vector<double, 5> weight5 = {0.23692689, 0.47862867, 0.56888889,
                                       0.47862867, 0.23692689};

struct Mxw {
  eg::MatrixXd k_matrix;
  eg::VectorXd x_list;
};

using Dd = std::function<double(double)>;
using Ddd = std::function<double(double, double)>;
class Gauss_IE_K_Matrix {
public:
  Gauss_IE_K_Matrix(int N) : m_N(N) {
    std::cout << __LINE__ << std::endl;
    switch (N) {
    case 5:
      gp = gaussp5;
      wg = weight5;
      break;
    case 7:
      gp = gaussp7;
      wg = weight7;
      break;
    case 9:
      gp = gaussp9;
      wg = weight9;
      break;
    case 11:
      gp = gaussp11;
      wg = weight11;
      break;
    default:
      std::cout << "invalid N = " << N << std::endl;
      exit(1);
    }
  }
  Mxw operator()(Ddd k, double a, double b, int n) {
    double h = (b - a) / n;
    eg::VectorXd x_list(m_N * n);
    eg::VectorXd w_list(m_N * n);
    eg::MatrixXd k_matrix(m_N * n, m_N * n);
    auto f = [&](int i) {
      double xi = i * h + a;
      double xf = (i + 1) * h + a;
      auto x = (xf + xi) / 2 + ((xf - xi) / 2 * gp).array();
      x_list.segment(i * m_N, m_N) = x;
      w_list.segment(i * m_N, m_N) = wg;
    };
    for (int i = 0; i < n; ++i) {
      f(i);
    }
    for (int i = 0; i < m_N * n; ++i) {
      for (int j = 0; j < m_N * n; ++j) {
        k_matrix(j, i) = w_list[i] * k(x_list[j], x_list[i]);
        // std::cout<< k_matrix(j,i) <<" ";
      }
    }

    return {k_matrix * h / 2, x_list};
  }

  int m_N = 0;
  eg::VectorXd gp;
  eg::VectorXd wg;
};

eg::VectorXd gauss_ie_sol(const Mxw &mxw, Dd h) {
  eg::VectorXd tmp(mxw.x_list.size());
  std::transform(mxw.x_list.cbegin(), mxw.x_list.cend(), tmp.begin(), h);

  eg::MatrixXd ik_inv =
      (eg::MatrixXd::Identity(mxw.x_list.size(), mxw.x_list.size()) -
       mxw.k_matrix)
          .inverse();
  return ik_inv * tmp;
}
